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Abstract 

We study the solution of the nonhnear BK evolution equation with the recently calcu- 
lated running coupling corrections [1,2]. Performing a numerical solution we confirm the 
earlier result of [3] (obtained by exploring several possible scales for the running coupling) 
that the high energy evolution with the running coupling leads to a universal scaling 
behavior for the dipole-nucleus scattering amplitude, which is independent of the initial 
conditions. It is important to stress that the running coupling corrections calculated 
recently significantly change the shape of the scaling function as compared to the fixed 
coupling case, in particular leading to a considerable increase in the anomalous dimension 
and to a slow-down of the evolution with rapidity. We then concentrate on elucidating 
the differences between the two recent calculations of the running coupling corrections. 
We explain that the difference is due to an extra contribution to the evolution kernel, 
referred to as the subtraction term, which arises when running coupling corrections are 
included. These subtraction terms were neglected in both recent calculations. We evalu- 
ate numerically the subtraction terms for both calculations, and demonstrate that when 
the subtraction terms are added back to the evolution kernels obtained in the two works 
the resulting dipole amplitudes agree with each other! We then use the complete running 
coupling kernel including the subtraction term to find the numerical solution of the re- 
sulting full non-linear evolution equation with the running coupling corrections. Again 
the scaling regime is recovered at very large rapidity with the scaling function unaltered 
by the subtraction term. 
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1 Introduction 



Recently our understanding of the linear Balitsky-Fadin-Kuraev-Lipatov (BFKL) [4, 5] and 
non-linear Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) [6-13] and 
Balitsky-Kovchegov (BK) [14-18] small-x evolution equations in the Color Glass Condensate 
[6-29] has been improved due to the completion of the calculations determining the scale of the 
running coupling in the evolution kernel in [1,2,30,31]. The calculations in [1,2] proceeded by 
including Nf corrections into the evolution kernel and by then completing Nf to the complete 
one-loop QCD beta-function by replacing Nf —Qt[(32- Calculation of the UsNj corrections 
is particularly easy in the s-channel light-cone perturbation theory formalism [32,33] used to 
derive the BK and JIMWLK equations: there as Nf corrections are solely due to chains of 
quark bubbles placed onto the s-channel gluon lines. 

The analytical results of [1,2] are not very concise and could not have been guessed without 
an explicit calculation. After finding «<, Nf corrections, the obtained contributions had to be 
divided into the running coupling part, which has a form of a running coupling correction 
to the leading-order (LO) JIMWLK or BK kernel, and into the "subtraction piece", which 
would bring in new structures into the kernel. Such separation had to be done both in [1] and 
in [2]. Unfortunately, there appears to be no unique way to perform this separation: it is not 
surprising, therefore, that it was done differently in both papers [1,2]. This resulted in two 
different running coupling terms, shown below in Eqs. fl35|) and fl36|) along with Eqs. ([7]) and 
([8]). Such a discrepancy has led to a misconception in the community that the calculations 
of [1] and of [2] disagree at some fundamental level. 

Indeed to compare the results of [1] and [2] one has to undo the separation into the running 
coupling and subtraction terms: combining both terms one should compare full kernels of the 
evolution equation obtained in [1,2]. There is another more physical reason to perform such 
comparison: in principle, there is no small parameter making the subtraction term smaller than 
the running coupling term and thus justifying neglecting the former compared to the latter. 
Even the labeling of one term as "running coupling" piece is somewhat misleading, since it may 
give an impression that the neglected subtraction term has no running coupling corrections in 
it. As was shown in [31] both terms actually contribute to the running coupling corrections to 
the BFKL equation (if one uses the separation of [2] to define the terms). 

In this paper we perform numerical analysis of the BK evolution equation with the Nf 
corrections resummed to all orders and with Nf completed to the QCD beta-function, Nf — > 
— 67r/52, with (32 given in Eq. fl2U]) . We first solve the BK equation keeping the running coupling 
term only, with the kernels given by Eqs. ([7j) and ([8]). Indeed the solutions we find this way 
are different from each other. We then evaluate the subtraction terms for both cases and show 
that inclusion of subtraction terms puts the results of [1] and [2] in perfect agreement with each 
other! We complete our analysis by solving the BK equation with the full kernel including both 
the running coupling and subtraction terms. 

This work is structured as follows. Section [2] begins with Sect. 12.11 in which we review 
the as Nf corrections to the dipole scattering amplitude evolution equation recently derived 
in [1,2] and the subtraction method employed in both works to separate the running coupling 
contributions from the subtraction terms. We discuss the scheme dependence of the running 
coupling terms introduced by this separation. We proceed in Sect. 12.21 by deriving the explicit 
expressions for the subtracted terms. The calculation is based on the results of [2]. Our 
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analytical results are summarized in Sect. 12.31 where we give the explicit final expression 
for the kernel of the subtraction term in Eq. fl39p . which, combined with Eq. fl38l) gives us the 
subtraction terms fl40l) and fl4ip for the subtractions performed in [1] and in [2] correspondingly. 

In Sect. [3] we explain the numerical method we use to solve the evolution equations. We 
also list the initial conditions used, along with the definition of the saturation scale employed. 
Throughout the paper we will avoid the important question of the Landau pole and the contri- 
bution of renormalons to small-x evolution. As we explain in Sect. [3], we will simply "freeze" 
the running coupling at a constant value in the infrared. For a detailed study of the renormalon 
effects in the non-linear evolution we refer the readers to [30]. 

Our numerical results are presented in Sect. |H By solving the evolution equations with the 
running coupling term only in Sect. 14.11 we show that the resulting dipole amplitude differs 
significantly from the fixed coupling case. We also observe that the amplitude obtained by 
solving the equation obtained in [2] is very close to the result of solving the BK evolution 
with a postulated parent-dipole running of the coupling constant. Both these amplitudes are 
quite different from the solution of the equation derived in [1], as one can see from Fig. |H 
In spite of that, all three evolution equations studied (the ones derived in [1], [2] and the 
parent-dipole running coupling model) give approximately identical scaling function for the 
dipole amplitude at high rapidity, as demonstrated in Fig. O in Sect. 14.21 It is worth noting 
that, as can be seen from Fig. [HI the anomalous dimension we extracted from our solution 
is 7 ~ 0.85, which is different from the fixed coupling anomalous dimension of 7 ~ 0.64. 
The former anomalous dimension also appears to disagree with the predictions of analytical 
approximations to the behavior of the dipole amplitude with running coupling from [34-38] . In 
Sect. 14.31 we numerically evaluate the subtraction terms for both [1] and [2] and show that their 
contributions are important, as shown in Fig. [71 However, subtraction terms decrease with 
increasing rapidity, such that at high enough rapidities their relative contribution becomes 
small (see Fig. [8]). In Fig. [9] we show that inclusion of subtraction terms makes the results of [1] 
and [2] agree with each other. Finally, the numerical solution of the full (all orders in as ^2) 
evolution equation including both the running coupling and subtraction terms is performed in 
Sect. 14. 4[ The results are shown in Fig. [TDi All the main features of the evolution with the 
running coupling are preserved in the full solution: the growth of the dipole amplitude and 
of the saturation scale with rapidity is slowed down (for the latter see Fig. [TT]) . The scaling 
function of Fig. [5] is unaltered by the subtraction term, as shown in Fig. [121 

We summarize and discuss our main conclusions in Sect. [3 

2 Scheme dependence 

2.1 Inclusion of running coupling corrections: general concepts 

The BK evolution equation for the dipole scattering matrix reads 



dY 



I 



d'^z K{xq,Xi,z) [S{xq,z; Y) S'(z, x^; Y) — S{xq,Xi; Y)] , 



(1) 



where 



K{xo,x^,z) 



(2) 
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is the kernel of the evolution. Here transverse two-dimensional vectors Xg and Xi denote the 
transverse coordinates of the quark and the anti-quark in the parent dipole, while z is the 
position of the gluon produced in one step of evolution [39-42]. We have introduced the 
notation r = Xjq — x^, = — z_, r_2 = ^~ ^Hi for the sizes of the parent and of the new 
(daughter) dipoles created by one step of the evolution. The notation r = |r| for all the 2- 
dimensional vectors will be also employed throughout the rest of the paper. Eq. ([1]) admits 
a clear physical interpretation: the original parent dipole, when boosted to higher rapidities, 
may emit a new gluon which, in the large- iVc limit, is equivalent to a quark-antiquark pair. 
Thus, the original dipole splits into two new dipoles sharing a common transverse coordinate: 
the transverse position of the emitted gluon, z_. The nonlinear term in the right hand side of 
Eq. ([1]) accounts for either one of the two new dipoles interacting with the target, along with 
the possibilities of only one dipole interacting or no interaction at all, while the subtracted 
linear term reflects virtual corrections. The kernel of the evolution is just the probability of 
one gluon emission calculated at leading logarithmic accuracy in ag^n-il/ xb), where xb is the 
fraction of momentum carried by the emitted gluon [39-42]. 

Under the eikonal approximation the dipole scattering matrix off a hadronic target at a 
fixed rapidity is given by the average over the hadron field configurations of Wilson lines V 
calculated along fixed transverse coordinates (those of the quark and of the antiquark). More 
specifically 

S{x^.x,-Y) = j^^{tv{V{x^)V\x^]). (3) 

Hence, the integrand of Eq. ([T]) can be regarded as a three point function in the sense that 
the gluon fields of the target are evaluated at three different transverse positions, those of the 
original quark and antiquark plus the one of the emitted gluon. 

However, the inclusion of higher order corrections to the evolution equation via all order 
resummation of Us Nj contributions as recently derived in [1, 2] brings in new physical chan- 
nels that modify the three point structure of the leading-log equation. The dipole structure 
generated under evolution by diagrams like the one depicted in Fig. [TJA. (for more detailed dis- 
cussion of the diagrammatic content of the high order corrections see [2]) is identical to the one 
previously discussed for the leading order equation, the only novelty being that the propagator 
of the emitted gluon is now dressed with quark loops, modifying the emission probability but 
leaving untouched the interaction terms. On the contrary, diagrams like the one in Fig. [TJ3 
in which a quark-antiquark pair (rather than a gluon) is added to the evolved wave function 
modify the interaction structure of the evolution equation. The evolution of the parent dipole 
scattering matrix driven by these kind of terms is proportional to the scattering matrix of the 
two newly created dipoles (the one formed by the original quark and the new antiquark and 
vice versa), ~ S{xq, Zi)S{z2^x_i). This term depends on four different transverse coordinates, 
i.e., it is a four point function and, therefore, its contribution to the evolution equation cannot 
be accounted for by a mere modification of the emission kernel of the leading order equation. 

To discuss in more detail the modifications introduced by the high order corrections, we 
find it useful to rewrite the evolution equation in the following, rather general way 

^-^^^^^ = :F[S{x„x,-Y)] (4) 

where ^ is a functional of the dipole scattering matrix which for the original derivation of the 
equation is given by the right hand side of Eq. ([1]). In general it can be decomposed into two 
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A B 

Figure 1: Schematic representation of the diagrams contributing to quark-NLO evolution, 
pieces 

j^[s] = n [s] -s[s] . (5) 

The first term, TZ, which we will call the 'running coupling' contribution, gathers all the 
higher order in Nj corrections to the evolution that can be recast in a functional form that 
looks identical to the leading order one but with a modified kernel, K, which includes all the 
terms setting the scale for the running coupling: 

T^[S{xq,x^;Y)] = j d^zk{xo,x^,z) [S{xo, z;Y) S{z, x^;Y) - S{xq,x^;Y)] . (6) 

The second term, S, henceforth referred to as the 'subtraction' contribution, encodes those 
contributions that depart from the three point structure of the leading-log equation. The 
explicit derivation and expressions for this term are presented in the next section. The relative 
minus sign between the two terms in Eq. ([5]) has been introduced for latter convenience. 

Importantly, the decomposition of JF into running coupling and subtraction contributions, 
although constrained by unitarity arguments, is not unique. Two different separation schemes 
have been proposed in [1,2]. They are both based on a similar strategy, sketched in Fig. O 
that can be summarized as follows. The newly created quark-antiquark pair added to the 
wave function in the diagrams Fig. [T]B is shrunk to a point, called the subtraction point, by 
integrating out one of the coordinates in the dipole-gg wave function, rendering the previously 
discussed four point nature of these contributions into a three point one. This integrated 
three point contribution is added to the running coupling contribution, whereas the original 
four point term minus its integrated version are assigned to the subtraction contribution. The 
divergence between the two approaches stems from the choice of the subtraction point. In 
the subtraction scheme proposed by Balitsky in [1] the subtraction point is chosen to be the 
transverse coordinate of either the quark, Z2, or the antiquark, Zi- The kernel for the running 
coupling functional, Eq. ([6]), obtained in this way is 



(7) 



On the other hand, in the subtraction procedure followed in [2] (which we will refer to as 
KW) the zero size quark-antiquark pair is fixed at the transverse coordinate of the gluon. 
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Figure 2: Schematic representation of the subtraction procedure. 



z = azi + (1 — a)z2, where a is the fraction of the gluon's longitudinal momentum carried by 
the quark, yielding the following expression for the kernel of the running coupling contribution: 



27r2 



'1 



as{rl) as{rl) r^-r^ , 



2 ^2 



rr r. 



1 ' 2 



where 



R^{r, r 1,7:2) 



ri r2 



(9) 



As we shall discuss later, the scheme dependence originated by the choice of the subtraction 
point is substantial and has an important effect in the solutions of the evolution equation when 
only the running contribution is taken into account. 

In our numerical study we will also consider the following ad hoc prescription for the kernel 
of the running coupling functional in which the scale for the running of the coupling is set to 
be the size of the parent dipole 

iV,a,(r2) r2 



27r2 



2 ^2 ■ 



rr r 



(10) 



1 ' 2 



This prescription is useful as a benchmark used to compare with previous numerical [3] and 
analytical works [34,35,43] where this ansatz was used. 



2.2 Derivation of the subtraction term 

We begin by considering the NLO contribution to the kernel of the JIMWLK and BK evolution 
equations with the s-channel gluon splitting into a quark- ant iquark pair, which then interacts 
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with the target, as shown on the left hand side of Fig. [31 The contribution of this diagram has 
been calculated in [2]. The resulting JIMWLK kernel is [2] 



da 



d'k d'k' d'q cPq' 



-ig- (z— Xg ■ [z—x-^ )~i(h~JSL ) '5.12 



X 



+ 



+ 



(27r)2 (27r)2 (27r)2 (27r)2 " 
1 

1 (1 - 2aYq ■ kk!_ ■ q' + q ■ q' k- y - q ■ y k- q' 



k^ + fa{l - a) y^ + g'2a(l - a 
2a(l-a)(l~2a) 



k:+fa{l-a) y^ + (f^a{l - a) 



k-q M. ■ q' 



12 



a 



^ + q^a{l 



a] 



k 



12 



q'^a[l 



a] 



The momentum labels in the above equation are explained on the left hand side of Fig. [31 
If k^ and k_2 the transverse momenta of the quark and of the antiquark in the produced 
pair as shown in Fig. [31 then the transverse momentum of the gluon is q = ki + k^- The 
other transverse momentum we use is = ^^^(l — a) — kj^a, where a is the fraction the of 
gluon's "plus" momentum carried by the quark, a = ki+/{ki+ + /c2+)- The prime over the 
transverse momentum denotes the momentum of the same particle in the complex conjugate 
amplitude. For instance q' is the momentum of the s-channel gluon in the complex conjugate 
amplitude. Finally, z-^ and Z2 denote the transverse coordinates of the quark and the antiquark. 
In Eq. (11 II) we use z^2 — ^1 ~^2 (^^^ transverse separation between the quark and the antiquark) 
and z = az^ + {l — a)z2 (the transverse coordinate of the gluon). 








^ 





Figure 3: A lowest order leading-A^j NLO correction which gives rise to the subtraction term 
is shown on the left. The same diagram with the gluon lines "dressed" by chains of fermion 
bubbles, as shown on the right, gives the full (resumming all powers of a^Nf) contribution to 
the subtraction term. Calculation of the subtraction term is pictured in Fig. [21 

To obtain the BK kernel from Eq. f[TT|) one should sum over all possible emissions of the 
gluon off the quark and antiquark lines in the incoming dipole both in the amplitude and in 
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the complex conjugate amplitude, which is accomplished by 

1 

m,n=0 

Below we will label the JIMWLK kernel by calligraphic letter /C and the corresponding BK 
kernel by K. 

The contribution of the kernel from Eq. f[T^ to the right hand side of the NLO version of 
Eq. is given by the following term 

"m j d^^id^^'2K^^°{xQ,x^;z^,Z2) S{xo,z^,Y) S{z2,x^,Y) (13) 

with the bare coupling. 

As shown in Fig. [2], at the NLO level, the subtraction term introduced in Eq. ([5]), is then 
defined by 

5nlo[5'] = "J y" d'^zid'^Z2K^^'^{xQ,x^;z^,Z2) 

X [5(xo,u7,r)5(^,xi,r)-5(xo,zi,r)5(z2,a,>^)], (i4) 

where w is the point of subtraction in the transverse coordinate space. In [1] it was chosen to 
be equal to the transverse coordinate of either the quark or the antiquark, 

w = z-^ or w = Z2, (15) 

as both choices lead to the same subtraction term ^^loI*^]- 

S^lolS] = I dh,dh2K^'^''{x^,x,;z,,Z2) 

X [Sixj^,z^,Y) S{z^,x^,Y) - S{xjQ,z^,Y) S{z2,x^,Y)]. (16) 
In [2] the subtraction point was chosen to be the transverse coordinate of the gluon z, 

w_ = z = a z^ + {1 — a) Z2. (17) 
This leads to the following subtraction term, which we denote ^^^fS*]: 

SlZiS] = I dh,dh2K^'^''{x„x,;z„Z2) 

X [S{a^,z,Y) S{z,x,,Y) - S{x^,z,,Y) S{z2,x,,Y)]. (18) 

Indeed the complete kernel in Eq. ([5]) is independent of the choice of w. However, since the 
subtraction term of Eq. (IT^ was neglected both in [1] and in [2] , different choices of w led to 
different expressions for the remaining running coupling part 7?.[S'], i.e., to different answers as 
far as investigations in [1] and in [2] were concerned. Different choice of w is the main source 
of the discrepancy of final answers of [1] and [2], though it does not imply any disagreement in 
the full expression ([5]). 
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Our goal in this Section is to evaluate Kf^^{xQ,Xi;zi,Z2) from Eq. (fTT!) including the 
running coupling corrections. The s-channel light-cone perturbation theory formalism makes 
such inclusion simple [2]: all we have to do is include infinite chains of quark bubbles on the 
gluon lines in the amplitude and in the complex conjugate amplitude, as depicted on the right 
hand side of Fig. [31 Performing calculations similar to those done in [2] one arrives at 



iNf I da 



d'k d''k' d'q d^q' 



(27r)2 (27r)2 (27r)2 (27r)2 
1 (1 — 2a)2g ■ k ■ q' -\- q ■ q' k- l£ — q ■ k- q' 



2n'2 



q^q 



k^ + g2a(l - a) k^^ + q'^a{l - a) 



2a(l-a)(l-2a) 



k-q ■ q' 



+ 



X 



k^ + q^a{l - a) k^^ + (f^a{l - a) 
4a2(i_a)2 



/2 



k^ + g2a(l - a) 



g^ + (fa{l 
1 



a 



qO. g 5/3 

1 + a^/?2 In =—z — 



1 + a^/?2 In - . 



(19) 



where 1C@ denotes the kernel with the running coupling corrections resummed to all orders. 
Just like in [2,31], here we will use MS renormalization scheme. Inclusion of fermion bubble 
chains generated two denominators at the end of Eq. (fT9ll . which is its only difference from 
Eq. ([n]). Here 



/?2 



llN^-2Nf 
12^? 



(20) 



Now we have to perform the transverse momentum integrals in Eq. (fT9l) . First we expand 
the denominators at the end of Eq. (JT9l) into a power series and rewrite Eq. ( fT9l) as 



^®(^0'^i;^i'^2) =4A^/ ^ (-aM/^2 



.n+m 



n,m=0 



dA" dX'- 



1 (1 - 2a)2g ■ ky ■ q' + q ■ q' k- y - q ■ y k- q' 




d'k d'k' d'q d^q' 
(2^(2^(2^(2^ 



^2 \ ^ / „i2\ ^' 



q2q/2 



+ 



k^ + g2a(l - a) 



2a(l-a)(l-2«) 



k'^ + q'^ail - 



a] 



k^ + q^a{l - a) k!^ + g'2a(l - 



a] 



k-q y ■ q' 



q^ 



q 



12 
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+ 



+ - a) 



a) 



(21) 



A=A'=0 



where we have defined = u^e^/^ to make the expressions more compact. 

Indeed we can not always expand the denominators of Eq. (fT9|) into a geometric series 
employed in Eq. (12T|) . but one has to remember that the summation of bubble chain diagrams 
shown on the right side of Fig. [3] gives one the geometric series. Hence the geometric series 
come first: later they are absorbed into the denominators shown in Eq. (JT9l) . which is an 
approximation not valid for all q and q'. Therefore, by keeping the geometric series in Eq. (!2Til 
we are not making any approximations. In general, in what follows we are not going to keep 
track of the issues of convergence of perturbation series. The contribution of renormalons to 
non-linear small-x evolution was thoroughly investigated in [30] and was found to be significant 
at low Q^. We refer the interested reader to [30] for more details on this issue. 

Using the following formulas 



d^k e - 



ik-z 



(27r)2 fc2 + g2 27r 



and 



d^k 



-ik-. 



k'^ + q^ 



(22) 



(23) 



we can now perform the k- and /c'-integrals in Eq. (12T1) . Integrating over the angles of q and q' 
as well yields 



ANt 



(27r) 



'^12 I I '2' ^1 



ra,m=0 

Aaaz 



rfA" d\' 



da dqqdq q { — 
J 





12 



, Z Xr\ I Z 



0) ±12 



U-^l) +^12 [Z-X, 



[z — 



X Ji{q\z — Xj^\) Ki{zi2 q ^faa) Ji{q' |^ — ^il) -^i(-2i2 Q.' \/oia) + 2aa{a — a) \/aa 
Ji{q \z - XqI) -^'1(2:12 q Vc^) Jo{q' \z - X]\) -^'0(^12 <?' ^faa) 



X 



^12 ' \^ 



Z12 \Z_ — Xj^\ 



1^' 



Z\2 \Z - ^1 



JoiQ \z - XqI) Kq{zi2 q \/oia) Ji{q' \z - x^\) Ki{zi2 q' Vaa) 



+ Aa^ Jo{q \ z-Xq\) Ko{zi2 q v^aa) Jo{q' \ z-Xi\) Kq{zi2 q' a/oq) 
We have defined 



A=A'=0 



(24) 



a = 1 — a 



(25) 
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for brevity. Now the integrals over q and q' can be carried out to give 



4 AT/ 
(2^ 



n,m=0 



da 



jJ^"^ oia 



A+A' 



r2(i + A) 



(1 + A) (1 + A') 



Z12 (tttt)^ 



X F 1 + A,2 + A:2: 



\Z 



F 1 + A',2 + A';2; 



X 



+ 



' ' — 2 

^12 ■ ~ i^o) n I \\ I ^ 1 \ o I \. o. U ~ i^ol 



^, -v2 



a a 2;: 



12 



+ 2 



« — a 



-12 



(1 + A)F 1 + A,2 + A;2; 



aa z 



12 



F 1 + A', 1 + A'; 1; 



' ' - 2 
a a 2:1 



^12 ■ — 1/ ri / 1 I \ 1 I \ . 1 . I— — ol 



^12 



F 1 + A,l + A;l; 



12 



(1 + A')F 1 + A',2 + A';2; 



a a 



+ 4-i^fl + A,l + A;l;- '- ) F ( 1 + AM + 1; 



^12 



' ' — 2 



(26) 



A=A'=0 



Unfortunately further simplification of the expression in Eq. fl26l) is impossible without approx- 
imations. The series resulting from summation over n and m are likely to be divergent due to 
renormalons. As we mentioned before, here we neglect the renormalon problem referring the 
reader to [30]. Similar to how it was done in [2] we are not going to attempt to resum the series 
exactly: instead we will calculate the next-to-leading order terms and assume that with a good 
accuracy they give us the scale(s) of the running coupling constant. This procedure is similar 
to the well-known prescription due to Brodsky, Lepage and Mackenzie [44]. 
Using the Taylor-expansions of hypergeometric functions 



F(l + A,2 + A;2;2) 



A 



1-2 l-z 



l + ln(l-z) + - ln(l-2) 

z 



o(A^ 



(27) 



and 



F(l + A,l + A;l;2) 
after some algebra we obtain 



A 

\-z l-z 



ln(l -2) +o(A2 



(28) 



^® (^0' ^11^7 



da 



« (^1 - ^0)^ + « (^2 - ^0)^] [« Ui - ^lY + « U2 - ^0)^] 42 



- Aaaz^2 - {z- Xq) z^2 " U - ^1) + ^12 U " ^0) ' U - ^1) 



1 - ttu /52 In 



2_ 

MS, 



oia, 



1 - "u /52 In 



2_ 

MS, 



oia, 
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+ 2 a a (a — a) z^2 i ^12 ' ii. ~ -^Lo) 

1 



1- a^P2 In 



1 



X 



X 



1 - tt/x /52 In 
1 - /?2 In 



2_ 

MS , 



1 

MS, 



12 



Rl{xQ)fi- 
l-a^P2 In 



+ oia. 



~^ ^12 ■ -ILi 
1 - /?2 In 



2_ 

MS, 



2_ 

MS, 



+ o a, 



I - a^[^2 In 



2_ 

MS, 



1 



Rl{x^)n- 



2_ 

MS, 



(29) 

In arriving at Eq. fl29l) we employed functions Rt{x) and Rl{x)j which have dimensions of 
transverse coordinates and are defined by 



In 



In 



4g-27-5/3 



[a {z-^ — xy + a{z2 — xy] ^ 



and 
In 



In 



aa ZI2 ^ 

{z — xy 

4g-27-5/3 



MS, 

« ~ 2;)^ + a U2 ~ 

^, ^2 



a a 2;: 



12 



(30) 



[tt (^1 - ^)2 + a (^2 - ^)2] 



-In 



0; (^1 — xY + a{z2 — x) 
a a zf2 



(31) 

The subscripts T and L stand for transverse and longitudinal gluon polarizations which give 
rise to the two different functions under the logarithm. 

Recombining the series in Eq. fl29l) into physical running couplings finally yields 



allCQ{xo,x^;z^,Z2) 



Nf 



da 



1 



47r'^ J [a {z_i — Xq)^ + a (^2 ~ i5.o)^] ["^ (^1 ~ i^Ii)^ + (^^2 ~ i^o)^] ^' 



21 .f^ 



4 a a ■ (^ - ^0) ^12 ■ - ^1) + ^12 U - ^0) ■ {z-x^] 



a 



Rxi^o) 



a. 



^ d2 



Rrp{X_l/ 



+ 2 a a (a — a) z 



12 



^12 ■ \^ ^0/ I d2 



Rrp{X_Qj 



a. 



+^12 ■ (^ ~ 3il) 



+ Aa'^ a'^ zfr, a 



12 I d2 



with the physical running coupling in the MS scheme given by 



as{l/R') 



a. 



1 + af,(32 In 



MS 



(32) 



(33) 



Eq. fl32|) is the contribution to the JIMWLK evolution kernel of the resummed diagram on 
the right hand side of Fig. O 
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2.3 Brief summary of analytical results 

Let us briefly summarize our analytical results. The non-linear small-x evolution equation with 
the running couphng corrections included reads 



dY 



n[S]-s [S] . 



(34) 



The first term on the right hand side of Eq. is referred to as the running couphng 
contribution. It was calculated independently in [1] and in [2]: the results of those calculations 
are given above in Eqs. (jT]) and (jSj) correspondingly, which have to be combined with Eq. ([6]) 
to obtain 



-Bal 



and 



^Kwj^j = d'zK'^'^{x^,x^,z)[S{xja,z-Y)S{z,x^]Y)- S{xs^,x^]Y)]. 



(35) 



(36) 



One notices immediately that 7^^^' [S] calculated in [1] is different from 7V^^ [S] calculated 
in [2] due to the difference in the kernels K^^^ and in Eqs. ([7]) and ([8]). However that 

does not imply disagreement between the calculations of [1] and [2]: after all, it is the full 
kernel on the right of Eq. flMl) . 7l[S] — S [S], that needs to be compared. To do that one has 
to calculate the second term on the right hand side of Eq. 

The second term on the right hand side of Eq. ( IMl) is referred to as the subtraction contri- 
bution. It is given by 

= j '^^^1 '^^^2 ^®(^o'^i;^i'^2) ['S'(^0':^)^) 'S'(w,Xi,F) - 5(^0,^1, F) 5'(^2,Xi,r)] 

(37) 

with the resummed BK kernel 

1 

K®{Xo,X^;Z^,Z2) = Cf (-l)™^"''C®fem,^n;^l,^2)- 



(38) 



m,n=0 



The resummed JIMWLK kernel /C@(x^, x„; z^, ^2) is given by Eq. fl32|) . along with Eqs. fl30l) 
and fl3T]) defining the scales of the running couplings. In the numerical solution below we will 
replace Nj —611^2 in its prefactor, obtaining 



a»^®(^o>^i;^i.^2) 



3^ 



da 



a {z^ - Xo)2 + a (^2 - Xo)2] [a {z^ - x^)'^ + a (^2 - Xj^y] zf^ 



4 a aii2 ■ (1 - ^0) I12 • U - a) + ^12 U - ^0) ■ U - ^1) 



a. 



+ 2 a a (a — a) z 



12 



+z 



±12 ' ULi) Ois 



^12 V— ^0/ '^s 



a. 



-Rt(^o) 



Rlix,] 



4 2:^2 



a. 



(39) 
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This substitution is the same as for all other factors of Nf. The same substitution was performed 
in [2] to calculate the running coupling term. In fact, as was shown in [31], the linear part of the 
subtraction term (calculated using the prescription of [2]) contributes to the running coupling 
corrections to the BFKL equation. Therefore, in that case, the factor of Nf in front of Eq. fl32l) 
is definitely a part of the beta-function. Hence the replacement Nf —Qti(32 is justified even 
in the subtraction term. Once again, in the numerical solution below we will use Eq. fl39|) along 
with Eq. ( l38l) in Eq. ( 137|) to calculate the subtraction term S[S]. 

Substituting w_ = (or, equivalently, w_ = Z2) in Eq. ( 1371) would yield the subtraction term 

^Bal^^j ^^2 J (fzi(fz2K(T){xQ,Xi;Z^,Z2) 

X [^(xo,^i,>^)^Ui,a,>^)-5(xo,Zi,F)5(^2,a^i,>^)] (40) 

which has to be subtracted from 7^^^' [S] calculated in [1] and given by Eq. to obtain the 
complete evolution equation resumming all orders of Og Nf in the kernel. 
Substituting w = z = azi + {l — a)z2in Eq. (137|) yields 

^Kw^^j =af^ J d'^zi(fz2K(T){xQ,Xi;z^,Z2) 

X [S{x,,z,Y) S{z,x,,Y) - S{x,,z,,Y) S{z2,x,,Y)] (41) 

which has to be subtracted from 7^^^ [S] calculated in [2] and given in Eq. fl5B]) again to obtain 
the complete evolution equation resumming all orders of Nf in the kernel. We checked 
explicitly by performing analytic calculations that the two evolution equations obtained this 
way agree at the NLO and NNLO. Below we will check the agreement of the two calculations 
to all orders by performing a numerical analysis of the solutions of these equations. 

The above discussion demonstrates that the separation of the evolution kernel into the 
running coupling and subtraction pieces, as done in Eq. flMl) . is somewhat artificial, and has 
no small parameter justifying one or another separation prescription. Therefore, the small-x 
evolution equation including all running coupling (or, more precisely, Nf) corrections should 
combine both terms in Eq. ([31]). Below we will solve such evolution equation numerically to 
obtain the full small-x evolution with the running coupling. 



3 Numerical setup and initial conditions 

In our numerical study we consider the translational invariant approximation in which the 
scattering matrix is independent of the impact parameter of the collision, i.e., S = S{r,Y). 
To solve the integro-differential equations, corresponding to the BK equation with running 
coupling we employ a second-order Runge-Kutta method with a step size in rapidity AY = 0.1. 
We discretize the variable \r\ into 800 points equally separated in logarithmic space between 
rmin = 10^® and r^ax = 50. Throughout this paper, the units of r will be GeV~^, and those 
of Qs will be GeV. All the integrals have been performed using improved adaptative Gaussian 
quadrature methods. The accuracy of this numerical method has been checked in [3] to be 
better than a 4% in all the r range. 
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We consider three different initial conditions for tlie dipole scattering amplitude, N{r, Y) 
1 — S{r,Y). The first one is taken from the McLerran-Venugopalan (MV) model [22,23]: 



N^^{r,Y = 0) = 1 - exp 



2/n'2 



(42) 



where a constant term has been added to the argument of the logarithm in the exponent in 
order to regularize it for large values of r. The other two initial conditions are given by 



N^^{r) = 1 - exp 



27 

sJ 



(43) 



with 7 = 0.6 and 7 = 0.8. These two last initial conditions will be referred hereinafter as AN06 
and AN08 respectively. The interest in this ansatz, reminiscent of the Golec-Biernat-Wusthoff 
model [45], is that the small-r behavior A^"^^ cx r'^"' corresponds to an anomalous dimension 
1 — 7 of the unintegrated gluon distribution at large transverse momentum. (AN labels initial 
conditions with anomalous dimension.) Our choices 7 = 0.6 and 7 = 0.8 can be motivated 
a posteriori by the observation that the anomalous dimension of the evolved BK solution for 
running coupling lies in between those two values and the one for the MV initial condition, 
7 ~ 1 (see Section 14.21) . Thus, the choice of distinct initial conditions allows us to better 
track the onset of the expected asymptotic universal behavior that is eventually reached at 
high energies and to study the influence of the pre-asymptotic, non-universal corrections to the 
solutions of the evolution equations. To completely determine our initial conditions, we set 
= 1 GeV at Y=0 in Eqs. (|42]) and (|43]) and put A = 0.2 GeV. Although Q'^ is normally 
identified with the saturation scale, our definition of the saturation scale through the rest of 
the paper will be purely pragmatical and given by the condition 

Ar(r = l/Q,(r),r) = «:, (44) 

with K = 0.5. We have checked that this choice of /t, albeit arbitrary, does not affect any of the 
major conclusions to be drawn in the rest of the paper. 

Finally, in order to avoid the Landau pole and to regularize the running coupling at large 
transverse sizes we stick to the following procedure: for small transverse distances r < r/^, with 
Tfr defined by ^^(l/rjj,) = 0.5, the running coupling is given by the one loop expression 

aJl/r^) = (45) 

with Nf = 3 and A = 0.2 GeV, whereas for larger sizes, r > rjr, we "freeze" the coupling at a 
fixed value = 0.5. A detailed study of the role of Landau pole in non-linear small-x evolution 
is given in [30]. 

4 Results 

In this Section, we discuss our numerical results and how they compare to previous numerical 
work and analytical estimates. 
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4.1 Running coupling 

Fig. m shows the solutions of the evolution equation when only the running coupling contribution 
is taken into account, i.e., neglecting the subtraction term in Eq. flMl) . for different initial 
conditions and for the three schemes considered in this work: Balitsky's, given by Eqs. ([7]) and 
f l35l) . KW, given by Eqs. ([8]) and fl36l) . and the the ad hoc parent dipole implementation of the 
running coupling, shown in Eq. ([TU]) . 




Figure 4: Solutions of the BK equation at rapidities Y=0, 5, 15 and 30 (curves are labeled from 
right to left) for the three running coupling schemes considered in this work: KW (solid line), 
Balitsky (dashed line) and parent dipole (dashed-dotted lines). The initial conditions are MV 
(top), AN08 (middle) and AN06 (bottom). 

As previously observed in [3,46], the most relevant effect of including running coupling 
corrections in the evolution equation is a considerable reduction in the speed of the evolution 
with respect to the fixed coupling case. This is a common feature of the different running 
coupling schemes studied here and of other phenomenological ones considered in the literature 
(a detailed comparison between the solutions for fixed couphng evolution and for parent dipole 
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running coupling can be found e.g. in [3]). This is not a surprising result, since a generic effect 
of the running of the coupling is to suppress the emission of small transverse size dipoles, which 
is the leading mechanism driving the evolution. 

However, despite this common feature of the running coupling solutions, significant dif- 
ferences are found between the solutions obtained under different schemes as we infer from 
Fig. m In particular, the evolution is much faster with the KW prescription than with that of 
Balitsky. Equivalently, the KW prescription yields a stronger growth of the saturation scale 
with rapidity /energy than Balitsky's. Moreover, the solutions obtained when the parent dipole 
prescription is used lay much closer to those obtained within the KW scheme than to the ones 
obtained when Balitsky's scheme is applied, contrary to what was suggested in [1]. As argued 
before, the differences observed in the solutions obtained using the two subtraction schemes 
are entirely due to neglecting the subtraction contribution and reflect the arbitrariness of the 
separation procedure. 

4.2 Geometric scaling 

It has been found in previous analytical [34, 36, 47] and numerical studies on the solutions 
of the BK equation at leading order [3,48-50] and for different heuristic implementations of 
next-to-leading order corrections [3,46], including the parent dipole prescription for the running 
coupling also considered in this work, that the solutions of the evolution equation at high enough 
rapidities are no longer a function of two separate variables r and Y, but rather they depend 
on a single scaling variable, r = rQs{Y). This feature of the evolution, commonly referred to 
as geometric scaling, is an exact property of the solutions for fixed coupling evolution due to 
the conformal invar iance of the leading-log kernel, and has become one of the key connections 
between the saturation based formalisms and the phenomenology of heavy ion collisions and 
deep inelastic scattering experiments [51-57]. 

It can be seen from Fig. [5] that the solutions of the BK equation with the running coupling 
terms discussed in the previous section also exhibit the property of scaling, in agreement with 
the analytical study carried out in [38], shown by the fact that the rescaled high rapidity 
solutions lay on a single curve which is independent of both the running coupling scheme and 
of the initial condition. The scaling behavior of the solution is observed in the whole r range 
studied in this work, including the saturation region, r > 1. The tiny deviations from a pure 
scaling behavior observed in Fig. [5] may be attributed to the fact that the full asymptotic 
behavior is reached at even larger rapidities {Y > 80, [3]) than those achieved by the numerical 
solution performed in this work. 

Remarkably, the scaling function for both KW and Balitsky's scheme coincides with the one 
obtained with the parent dipole prescription, up to the above mentioned scaling violations. It 
has been observed in [3,46,48] that the scaling function differs significantly in the fixed and 
running coupling cases. Following that work, and to make a more quantitative study of the 
scaling property, we fitted our solutions to the functional form [34] 

/(r) =aT^^ {Inr^ + b) , (46) 

with a, b and 7 free parameters, within a fixed window below the saturation region, r G 
[10-^0.1]. Noticeably, at large enough rapidities the whole fitting window lays within the 
geometric scaling window proposed in [36]: {A/Qs{Y)) < r < 1, where A is some initial scale. 
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Figure 5: Solutions of the BK equation at rapidities Y=0 and 40 for KW (solid line), Balitsky 
(dashed line) and parent dipole (dashed-dotted lines) schemes plotted versus the scaling variable 
r = rQs(y). The asymptotic solution obtained with fixed coupling «>, = 0.2 at F = 40 in [3] 
is shown (black dashed-dotted line) for comparison. The initial conditions are MV (left) and 
AN06 (right). 



The value of 7 extracted from the fits at rapidity F = 40 lays in between 7 ~ 0.8 and 7 ~ 0.9. 
This conclusion holds for the three initial conditions used here: the anomalous dimension 
seems to converge to some intermediate value, in agreement with the value found in [3], for 
asymptotic running coupling solutions (7 ~ 0.85 at F = 70). This result for anomalous 
dimension is very far away from the value obtained in [3] for fixed coupling solutions (7 ~ 0.64 
at y = 70) and from the predicted anomalous dimension for both running and fixed coupling 
solutions from analytical studies of the equation based on saddle point techniques [34-38], 
7c = x(7c)/x'(7c) = 0.6275, where x is the leading-log BFKL kernel. 

It might be argued that the numerical value of the anomalous dimension extracted from our 
fits is conditioned by the choice of the fitting function and by the fitting interval. Actually, it 
was shown in [3] that the solutions of the evolution could be well fitted by other functional forms, 
including the double-leading-log solution of BFKL, within a similar fitting region to the one 
considered in this work. On the other hand, several phenomenological parameterizations of the 
solution of the evolution have been proposed in [54-57] and have successfully confronted HERA 
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Figure 6: Asymptotic solutions (Y=40) of the evolution equation for running coupling (solid 
line) and fixed coupling with = 0.2 (dashed line). A fit to a power-law function ar^'^ in the 
region r G [10~®, 10~^] yields 7 ~ 0.85 for the running coupling solution and 7 ~ 0.6 for the 
fixed couphng one. 



and RHIC experimental data. There, the dipole scattering amplitude at arbitrary rapidity is 
assumed to be given by a functional form analogous to our ansatz for the initial condition 
Eq. (H3ll . but allowing for geometric scaling violations by replacing 7 —>■ 7(r, F). The value 
of the anomalous dimension at r = 1/Qs and/or for y — > 00 is fixed to be the BFKL saddle 
point, 7c ~ 0.63 (the saddle point value considered in [57] is slightly different, 7 ~ 0.53), while 
the value 7 = 1 is recovered in the limit r — > 00 at any finite rapidity. The success of these 
phenomenological works supports the claim that the anomalous dimension of the solution is 
given by the BFKL saddle point, in agreement with the above mentioned analytical predictions. 
However, the relevant values of momenta probed at current phenomenological applications are 
very distinct from the fitting region considered here. For example, the inclusive structure 
function measured in HERA is fitted in [54,56] within the region 0.045 GeV^ < < 45 GeV^, 
whereas charged hadron pt spectra in dAu collisions is well reproduced by [55-57] in the region 
1 GeV < < 4.5 GeV. Note that, for both sets of data, the measured regions overlap with 
the deeply saturated domain of the solution. On the contrary, our fitting region 10^^ < r < 
1 corresponds to values of momenta ~ 10 Qs{Y) < pt < 10^ Qs{Y) (always well above the 
saturation scale), with Qs{Y = 40) ~ 500^1000 GeV for the different running coupling schemes 
considered and, therefore, has no overlap with the kinematic regions measured experimentally, 
since we scrutinize a momentum region strongly shifted to the ultraviolet compared to currently 
available data. Moreover, it should be noticed that the rapidity interval covered by both sets 
of experimental data is AY < 4 in both cases, while we study the solutions of the evolution 
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Figure 7: Subtraction contribution calculated in the KW scheme (triangles) and in Balitsky's 
(stars). The trial functions correspond to the solutions of the evolution under Balitsky run- 
ning coupling scheme at rapidities Y = 0,30 for MV (left), AN08 (center) and AN06 initial 
conditions. 



at asymptotic rapidities, Y ~ 40. We have checked that shifting our fitting region to larger 
values of r (smaller momentum) would bring the value of 7 extracted from our fits closer to the 
saddle point BFKL one, since the transition from the ultraviolet region to the deeply saturated 
domain of the scaling solution is realized by a locally less steeper function (see Figs. ([5]) and 
f ll2l) ). Therefore, there is no contradiction at all between the success of the phenomenological 
parameterizations of the solutions and the results reported here. 

With the above clarifications we reach the following conclusion: the asymptotic scaling 
solutions corresponding to fixed and running coupling evolution are intrinsically different in 
the whole r-range. This is emphasized in Fig. [6], where we represent the scaling solutions in 
a log scale for r < 1. It is clear that the tail of the distribution falls off with decreasing r 
much steeper for the running coupling solution than for the fixed coupling one. A fit to a 
pure power-law function, / = ar^'*', in the region r G [10^^,10^^] yields 7 ~ 0.85 for the 
running coupling and 7 ~ 0.61 for the fixed coupling solution. The differences between fixed 
and running coupling solutions at r > 1 are evident from Fig. [51 This is a puzzling result that 
remains to be understood from purely analytical methods. 
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Figure 8: Ratio of the subtraction over the running terms, V{r,Y) = S[N{r,Y)]/TZ[N{r,Y)], 
calculated in both KW (triangles) and Balitsky (stars) schemes for MV (left), AN08 (middle) 
and AN06 (right) initial conditions at rapidities Y=0 (top) and Y=30 (bottom). 

4.3 Subtraction Term 

Before attempting to solve the complete evolution equation, and in order to gain insight in the 
nature and structure of the subtraction contribution, we first evaluate the subtraction functional 
for both Balitsky, Eq. (HOj) . and KW, Eq. (l4T|l . schemes using a set of trial functions for S which 
we choose to consist of the solutions of the evolution equation with the running coupling in 
Balitsky's scheme at different rapidities and of the three initial conditions considered above in 
this work. 

Two main remarks can be made about our results, shown in Fig. [T) 

i) For all the trial functions considered in this work, the subtraction contribution is much 
larger in the KW scheme than in Balitsky's. A plausible explanation for this is that Balitsky's 
subtraction contribution, Eq. fHOj) . when expanded in terms of dipole scattering amplitudes, 
N = 1 — S, reduces to a sum of non-linear terms, since all the linear terms in the expansion 
cancel each other due to the Zi ^ z_2 symmetry of the kernel, whereas in the KW case no such 
cancellation happens and the subtraction contribution, Eq. (HTI) . also includes linear terms, 
which are dominant over the non-linear ones in the non-saturated domain where <^ 1. 

ii) The subtraction contribution S has the same sign as the running coupling contribution 
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TZ in the whole r range which, together with the relative minus sign assigned to the subtraction 
term in Eq. implies that the proper inclusion of the subtraction term reduces the value 

of the functional that governs the evolution, JF. In other words: the subtraction contribution 
tends to systematically slow down the evolution, as we shall explicitly confirm in the next 
subsection. 

To better quantify the size of the subtraction contribution, we plot the ratio 'D{r, Y) = 
S[N{r, Y)]/TZ[N{r, Y)] in Fig. [HI At Y = 0, the relative weight of the subtraction contribution 
with respect to the running one within the KW scheme and for a MV initial condition goes 
from a P ~ 0.4 at small rtoP~latr~l. The same ratio for the Balitsky scheme takes 
significantly smaller values: it goes from P ~ 0.1 at small r to ~ 0.4 for r ~ 1. As the 
evolved solutions get closer to the scaling function, i.e. for larger rapidities, the r dependence of 
the ratio becomes flatter and its overall normalization goes down to an approximately constant 
value V ~ 0.15 for the KW scheme and T> ~ 0.025 for that of Balitsky. This behavior remains 
unaltered when going from rapidity F = 20 to F = 30, which suggests that the ratio may 
saturate to a fixed value in the asymptotic region. 




Figure 9: Total kernel T = TZ — S calculated under Balitsky's scheme, Eqs ([7]) and ( HOl) . (solid 
line) and under the KW scheme, Eqs ([H]) and ( HTl) . (dashed line). The overlap of the two lines 
shows the agreement between the two calculations. Triangles stand for the running coupling 
term calculated in the KW approach, TZ^^ , while stars stand for the running coupling term 
under Balitsky's scheme, 7^^*^^ The trial functions N{r, Y) correspond to the solution of the 
evolution with only running coupling under Balitsky's scheme at Y=0 (left) and Y=30 (right) 
for a MV initial condition. 
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Figure 10: Solutions of the complete (all orders in ag P2) evolution equation given in Eq. fl34l) 
(solid lines), and of the equation with Balitsky's (dashed lines) and KW's (dashed-dotted) 
running coupling schemes at rapidities y = 0, 5 and 10. Left plot uses MV initial condition. 
The right plot employs the initial condition given by the dipole amplitude at rapidity Y = 35 
evolved using Balitsky's running coupling scheme and with r-dependence rescaled down such 
that Qs = Q's = l GeV. 



Finally, we have checked that combining the subtraction and running coupling contributions 
for both schemes adds up to the same result. This is shown in Fig. [9l where we plot the value 
of the total functional T = IZ — S calculated under the KW scheme (Eqs. (ED and ( l36l) for 
the running coupling term, TZ, and Eq. (HT!) for the subtraction term, S) and under Balitsky's 
scheme (Eqs. ([7]) and (!35|l for the running coupling term and Eq. (HOl) for the subtraction term). 
The two results coincide within the estimation of the numerical accuracy previously discussed. 
The agreement between the two results is better in the small-r region, where the two curves lay 
almost on top of each other. In the saturation region, r > 1, the agreement is slightly worse, 
although the differences between the values of JF calculated in both schemes is still much less 
than the differences between the running coupling terms themselves. This slight remaining 
disagreement between the Balitsky's and KW prescriptions may also be due to inaccuracies in 
a Fourier transform of a geometric series performed in arriving at Eq. (l39|) . This result serves 
as a cross-check of our numerical method and as an additional confirmation of the agreement 
of the independent calculations derived in [1,2]. 
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Figure 11: Saturation scale corresponding to the solutions plotted in Fig. [TOl 



4.4 Complete running coupling BK equation 

In this section we calculate the solutions of the complete evolution equation, Eq. including 
both the running and subtraction terms obtained by the all-orders Nj resummation and by 
the Nf — * —67^(32 replacement. Since the numerical evaluation of the subtraction contribution 
at each point of the grid and each step of the evolution would require an exceedingly large 
amount of CPU time consumption, the strategy followed to include it in the evolution equation 
consists of calculating such contribution only in a small set of grid points at each step of the 
evolution, which we fixed at n = 16, between the points ri and r2, which are determined at each 
step of the evolution by the conditions N{Y,ri) = 10~^, and N{Y,r2) = 0.99, and then using 
power-law interpolation and extrapolation to the other points of the grid. Both the running 
and subtracted terms are calculated within Balitsky scheme. This procedure is motivated by 
the fact that, as discussed in the previous section, the subtraction contribution can be regarded 
as a small perturbation with respect to the running coupling term within Balitsky's scheme and 
by the fact that it is a rather smooth function that can be well fitted by a power-law function in 
most of the r-range. The accuracy of this procedure has been checked by doubling the number 
of points at which the subtraction contribution is calculated at each step of the evolution, i.e. 
by setting n=32. At Y=2, the differences between the solutions obtained with the two above 
mentioned choices for n were less than a 8% in the tail of the solution, r < ri, and less than a 
3% for r > ri. 

The results of the evolution calculated in this way and using MV and rescaled asymptotic 
running coupling solution (Y=35) as initial conditions are plotted in Fig. (TDl They confirm 
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10"^ 10"* 1 ^ 10 



Figure 12: Rescaled solutions given by the complete /52-evolution equation (solid line) and for 
KW (dashed-dotted line) and Balitsky's (dashed line) running coupling schemes at F = 10. The 
initial condition corresponds to the dipole amplitude at rapidity y = 35 evolved using Balitsky's 
running coupling scheme and with r-dependence rescaled down such that Qs = Q'^ = 1 GeV. 



the expectations raised in the previous Subsection: the inclusion of the subtraction terms 
considerably slows down the evolution with respect to the sole consideration of the running 
coupling contributions. Moreover, the reduction in the speed of the wave front is much larger 
for the KW scheme than for that of Balitsky one for both initial conditions. However, the closer 
the initial condition is to the asymptotic running coupling scaling function, the smaller are the 
effects of the subtraction contribution. These features can be better quantified by inspecting 
the rapidity dependence of the saturation scale generated by the evolution, plotted in Fig. [TTl 
At rapidity F = 10 the ratio of the saturation scale Qs yielded by the KW scheme to Qs given 
by the complete /32-evolution equation is a factor of ~ 2.5 for the MV initial condition and 
a factor of ~ 2.1 for the asymptotic running coupling initial condition. At the same rapidity, 
the ratio of the saturation scale obtained under Balitsky's scheme to Qs corresponding to the 
complete /52-evolution is ~ 1.25 for the MV initial condition and ~ 1.15 for the scaling 
function initial condition. Thus, in spite of the smallness of the ratio of the subtraction terms 
to the running coupling contributions at high rapidity, which is ~ 0.025 for Balitsky's and 
~ 0.15 for KW scheme at F = 30 (see bottom plots in Fig. [H]), the proper inclusion of the 
subtraction term results in fairly sizable effects in the solutions of the evolution equation. 

Finally, we notice that the scaling behavior of the solution is not affected by the subtraction 
term. This is seen in figure Fig. [121 where we evolve starting from an initial condition already 
close to the running coupling scaling function and plot the solutions of the evolution equation 
obtained with just running coupling terms (see Section 14.11) and the solution of the complete 
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ds /32-evolution at rapidity Y = 10. It is clear that, within the numerical accuracy, no departure 
from the scaling behavior is observed. Therefore the main effect of a proper consideration of the 
subtraction term is the one of reducing the speed of the evolution. It does not violate or modify 
the geometric scaling property of the solutions established in Section 14.21 In our understanding 
geometric scaling appears to persist when the running coupling effects are included because, 
at high enough rapidity Qs{Y) ^ A, such that the new (from the LO standpoint) momentum 
scale A introduced by the running coupling can be safely neglected. Hence the dynamics is 
again characterized by a single momentum scale Qs{Y). At the same time running coupling 
does modify the evolution kernel, leading to a different shape of the scaling function. 

5 Conclusions 

In this paper we have taken into account all corrections to the kernels of the non-linear JIMWLK 
and BK evolution equations containing powers of agNj. We reiterated the fact that the sep- 
aration of the resulting kernel resumming all powers of Nf into the running coupling and 
subtraction parts, as done in the previous calculations of [1,2], is not justified parametrically. 
We have then performed numerical analysis with the following conclusions. 

• First we solved the evolution equations derived in [1] and [2] keeping only the running 
coupling part or the evolution kernel and neglecting the subtraction term. Comparing to 
the results for fixed coupling obtained in [3] we confirmed the conclusion reached in [3] 
that the growth with rapidity is substantially reduced when running coupling corrections 
are included. The results for three different initial conditions are shown in Fig. |H We 
observe that the solution of the equation derived in [2] differs significantly from that 
derived in [1], but agrees (with good numerical accuracy) with the solution of the BK 
evolution equation with the coupling running at the parent dipole size. (The latter is just 
a model of the running coupling not resulting from any calculations, which we plot for 
illustrative purposes.) We also observe that at sufficiently high rapidity both equations 
from [1] and from [2] give us the same scaling function for the dipole amplitude N{r, Y) 
as a function of rQs{Y), which is also in agreement with the scaling function given by 
the parent dipole running, as shown in Fig. [5l The fact that the scaling is preserved when 
the running coupling corrections are included was previously established in [3], though 
for models of running coupling only. The shape of the scaling function is very different 
from that obtained from the fixed coupling evolution equations. In particular, we found 
that for dipole sizes below 0.1 /Q^ the anomalous dimension of the scaling function in the 
running coupling case becomes 7 ~ 0.85 (see Fig. [6]). This is different from the result of 
several analytical estimates [34-38], which expect the anomalous dimension not to change 
when running coupling corrections are included and to remain at its fixed coupling value 
of 7 ^ 0.63. 

• We have then evaluated the subtraction term for both calculations performed in [1] and 
[2]. We demonstrated that subtracting the subtraction terms from the running coupling 
terms makes the full answer agree for both calculations of [1] and [2] , as shown in Fig. M 
for the right hand side of the evolution equation. It turns out that the subtraction 
term 4S^'^^[5'], which has to be subtracted from the result of [1], is systematically smaller 



25 



than iS^^[S'], to be subtracted from the result of [2], over the whole rapidity range 
studied here. This implies that the result of [1] should have a smaller correction than 
the result of [2] and is thus closer to the full answer. The subtraction terms iS^^'fS*] 
and [5] are plotted in Fig. [7] as functions of the dipole size r for different values of 
rapidity. Their relative contributions to the evolution kernel are shown in Fig. [8|, where 
we plotted the subtraction functional divided by the running coupling functional. From 
those figures we conclude that both the magnitude of these extra terms and their relative 
contribution to the evolution kernel decrease with increasing rapidity. Hence, while at 
"moderate" rapidities (the ones closer to realistic experimental values) the subtraction 
term is important for both calculations [1,2], it becomes increasingly less important at 
asymptotically large rapidities. The physics is easy to understand: the subtraction terms 
are o{al), while the running coupling part of the kernel is o{as). Hence, if we suppose 
that the effective value of the coupling is given by its magnitude at the saturation scale 
Qs{Y), then, as rapidity increases, the couphng would decrease, making the subtraction 
term much smaller than the running coupling term. Indeed, while at asymptotically high 
rapidities the assumption of [1,2] that the subtraction term could be neglected is justified, 
making the results of [1] and [2] agree with each other, for rapidities relevant to modern 
days experiments the subtraction term is numerically important. 

• With the last conclusion in mind we continued by numerically solving the full evolution 
equation resumming all powers of Nf in the evolution kernel, which now would combine 
both the running coupling and the subtraction terms. The five- dimensional integral in 
the subtraction term (l37j) made obtaining this solution rather difficult. The outcome of 
the calculation is shown in Fig. [TOl All the main conclusions stated above were again 
confirmed by the solution of the full equation. At asymptotically high rapidity scaling 
regime is recovered, as can be seen from Fig. [121 As the subtraction term is less important 
in that regime, the scaling function appears to be the same as in the case of having only 
the running coupling term in the kernel. The anomalous dimension again turns out 
to be 7 ~ 0.85, in disagreement with the analytical expectations of [34-38]. However, 
the scaling of the saturation scale with rapidity appears to be in agreement with the 
expectations of analytical work of [34,35,38], as shown in Fig. [TTl 

We conclude by observing that the knowledge of the non-linear small-x evolution equation 
with all the running coupling corrections included brings us to an unprecedented level of pre- 
cision allowing for a much more detailed comparison with experiments than was ever possible 
before. 
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